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Abstract. We propose a finite elements algorithm to solve a fourth order partial dif- 
ferential equation governing the propagation of time-harmonic bending waves in thin 
elastic plates. Specially designed perfectly matched layers are implemented to deal 
with the infinite extent of the plates. These are deduced from a geometric transform in 
the biharmonic equation. To numerically illustrate the power of elastodynamic trans- 
formations, we analyse the elastic response of an elliptic invisibility cloak surround- 
ing a clamped obstacle in the presence of a cylindrical excitation i.e. a concentrated 
point force. Elliptic cloaking for flexural waves involves a density and an orthotropic 
Young's modulus which depend on the radial and azimuthal positions, as deduced 
from a coordinates transformation for circular cloaks in the spirit of Pendry et al. [Sci- 
ence 312, 1780 (2006)], but with a further stretch of a coordinate axis. We find that 
a wave radiated by a concentrated point force located a couple of wavelengths away 
from the cloak is almost unperturbed in magnitude and in phase. However, when the 
point force lies within the coating, it seems to radiate from a shifted location. Finally, 
we emphasize the versatility of transformation elastodynamics with the design of an 
elliptic cloak which rotates the polarization of a flexural wave within its core. 
AMS subject classifications: 52B10, 65D18, 68U05, 68U07 

Key words: Finite Elements; High-Order Differential Operator; Perfectly Matched Layers; Trans- 
formation Elastodynamics; Cloaking 



1 Introduction 

In 2006, Pendry et al. have shown that by surrounding a finite size object with a coating 
consisting of a metamaterial, we can render this system transparent to electromagnetic 
radiation [IJ. The cornerstone of this work is the geometric transformation 
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where r' , 6' and z' are radially contracted cylindrical coordinates and (x,y,z) is the Carte- 
sian basis. This transformation maps the disk Dr^ \ {0} = {'':0<r<R2} onto an annulus 
Dr2\^Ri = • ^1 ^ ^ ^2}- In other words, if a source located in R^\Dr2 radiates in 
vacuum, the electromagnetic field cannot reach the disk Dr^ and therefore this region 
is a shelter for any object. In layman terms, a point in the space is transformed into a 
disc leading to the creation of a hole in the physical space where waves cannot penetrate 
from outside but are guided around this area in a way that anything we put inside this 
disc would be invisible to exterior observers, in the context of electromagnetism [Ij-El, or 
neutral, in the context of elastodynamics ||6HT2|. 

From a geometric point of view, the change of coordinates means that in the annulus 
we should work in a stretched space with associated metric tensor 

T = J^,J„Vdet(J„0, (1-2) 

d(v z) 

with Jrr' = ^(j/q! z') Jacobian of the transformation from free to contracted cylindrical 
coordinates, its transpose and det(Jr]./) its determinant. 

In terms of material parameters, the only thing to do in the annulus is to replace the 
material (homogeneous and isotropic) by an equivalent one that is inhomogeneous (its 
characteristics are no longer piecewise constant but merely depend on r', 6', z coordi- 
nates) and anisotropic ones (tensorial nature) whose properties are given by ||3l [T3HT5| 

^=eT-\ and ^^=^^-^. (1.3) 

Variants of this formula appear in many instances in the existing literature, such as in 
the book on the geometry of electromagnetism by Post, back in 1962 ||T6|. However, it 
seems that Pendry and Ward were the first ones to apply it to simplify the computational 
analysis of certain types of photonic materials in 1996 [171 . 

The diagonalised form of the tensors £^ = Diag(£r,£0, £3) and = Diag(/^r,//0,//3) was 
given in m 

r-Ri r f R2 \'r-Ri 

^r = Hr = ^,ee = ^e = j^^,es = ^,= ^-^-^J (1.4) 

Mimicking the heterogeneous and anisotropic nature of the tensors of permittivity 
and permeability given by eq. H1.4I) became possible only when Pendry suggested the 
use of newly discovered metamaterials which are composite structures manufactured for 
their exotic properties enabling one to control the electromagnetic field. A team lead 
by Pendry and Smith implemented this idea using a metamaterial consisting of concen- 
tric layers of Split Ring Resonators (SRR), which made a copper cylinder invisible to an 
incident plane wave at 8.5 GHz as predicted by the numerical simulations [2J. Inde- 
pendently, Leonhardt studied conformal invisibility by solving the Schrodinger equation 
which is valid in the geometric optics limit |fl8| . 
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A different path to invisibility was followed by McPhedran et al. flS] who proposed 
to cloak a countable set of line sources using anomalous resonance when it lies in the 
close neighbourhood of a cylindrical coating filled with a negative index material which 
is nothing but a cylindrical version of the celebrated perfect lens of Sir John Pendry ||20| . 

An impedance matching route to invisibility was further developed by Alu and En- 
gheta [21 J. Although very promising, their proposal relies on a specific knowledge of 
the material properties of the object being concealed. Recently, these authors gave the 
experimental proof of their idea p2| . 

Other routes to invisibility are based on homogenization; A team led by Shalaev HI 
has shown the possibility to make an object nearly invisible in TE polarization. A reduced 
set of material parameters was introduced to relax the constraint on the permeability, thus 
leading to an impedance mismatch with vacuum which was shown to preserve the cloak 
effectiveness to a good extent. 

Farhat et al. f231 analyzed cloaking of transverse electric (TE) fields through homoge- 
nization of radially symmetric metallic structures. The cloak consists of concentric layers 
cut into a large number of small infinitely conducting sectors which is equivalent to a 
highly anisotropic permittivity. This structure was shown to work for different wave- 
lengths provided they are ten times larger than the outermost sectors. 

Soon after. Cummer and Schurig 1 6 [ analyzed the 2D acoustic cloaking for pressure 
waves in a transversely anisotropic fluid by exploiting the analogy with TE electromag- 
netic waves. Torrent and Sanchez-Dehesa \ 7\ subsequently investigated this cloaking for 
concentric layers of solid lattices behaving as anisotropic fluids in the homogenization 
limit. Using a similar approach, Farhat et al. |8| independently demonstrated cloaking of 
surface liquid waves using a micro-structured metallic cloak which was experimentally 
validated at 10 Hz. Quite remarkably, Chen and Wu [9J and Cummer et al. fM] noticed 
that a 3D acoustic cloaking for pressure waves in a fluid can be envisaged since the wave 
equation retains its form under geometric changes. 

However, when one moves to the domain of elastodynamics, Milton, Briane and 
Willis 1 10] have shown that the governing equations are not invariant under coordinate 
transformations and consequently that if cloaking exists for such type of waves, it would 
be of a different nature than its acoustic and electromagnetic counterparts. Indeed, the 
equation of propagation of elastic waves in solids (Navier equation) with a time harmonic 
dependence can be written in the weak form: 

V-C:Vu+po<^^u = 0, (1.5) 

where u{x,y,z,t) = u{x,y,z)e~'^\ is the three-component vector displacement field, po is 
the scalar density of the elastic medium, C is the rank 4 elasticity tensor, co is the wave 
angular frequency, and t is the time. 

Moreover, Milton, Briane and Willis noticed IfTOll that by introducing the geometric trans- 
form X— >x' such that u'(x')=A-^u(x) with Aij=dx'j/dXj, thereby enforcing the symmetry 
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of the elasticity tensor, equation (|1.7|) takes the form: 



V • (C + S') : Vxi'+p'M = D' : V'u' , 



(1.6) 



which importantly preserves the symmetry of the new elasticity tensor C'+S'. However, 
this transformed equation contains S' and D' which are two rank 3 (symmetric) tensors 
such that Dpg,. = S^,.p and p'^^ is a rank 2 tensor whose expressions can be found in [lOJ. 

Nevertheless, Brun, Guenneau and Movchan recently discovered that if one does not 
constrain the symmetry of the elasticity tensor, equation il.7i takes the form for the case 
of fully coupled in-plane shear and pressure waves: 



where C is a rank four tensor with only the main symmetries and p' a scalar quantity, 
both of them spatially varying, see 

Farhat, Guenneau, Enoch and Movchan further considered the case of biharmonic 
bending waves propagating in thin plates and have shown that by considering a radi- 
ally dependent isotropic mass density and a radially dependent and orthotropic flexural 
rigidity, it was possible to extend cloaking by refraction to the case of out-of-plane elastic 
waves [25 J. 

In this paper, we give a numerical analysis of the biharmonic problem based on the 
Finite Elements Method (FEM) which is an approximation technique for the partial differ- 
ential equations (PDE) and require a variational formulation of the problem to be stud- 
ied (expressing expressions in a weak form). Indeed, the COMSOL formulation of the 
fourth order biharmonic equation is given, supplied with appropriate boundary condi- 
tions (clamped or stress-free) and Perfectly Matched Layers (PMLs). In the last section, 
elliptic cloaking is numerically performed and the rotator and mirage effects established 
for these waves confirming that they behave in a way similar to electromagnetic and 
acoustic waves. 

2 Finite Element modelling for the biharmonic equation 

The equations for bending of plates are well known and can be found in many textbooks, 
such as those of Timoshenko or Graff [26J- [27]. The displacement is W{x,y,t) in the z- 
direction. We choose to work in cylindrical coordinates with a time harmonic dependence 
i.e. u{x,y,t) = (0,0,W(x,t/)e^'P'), where t is the time and p is the angular frequency. 

The wavelength A is supposed to be large enough compared to the thickness of the 
plate h and small compared to its in-plane dimensions h <^ \ <^ L. In this case we can 
adopt the hypothesis of the theory of Von-Karman. 

To derive the equation of motion of the flexural out-of-plane waves, we can use a 
variational point of view pS] . 



V-C: Vu+pVu = 0, 



(1.7) 
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2.1 Weak form associated with the fourth-order equation 

We now consider the following fourth order elliptic equation 

flV-(bV(flV-(bVW)))=/ in Q, (2.1) 

where a is a scalar, and b is a matrix describing the material of the domain, and / a source 
term. This equation is supplied with 'Dirichlet' (or clamped) boundary conditions m = 
and n- VW = on 3Q (we note that the second boundary condition is needed since we 
are in presence of a fourth-order partial differential equation). 

By multiplying Equation (|2.ip with a test function W and integrating on the domain O 

/ dTr(V-(bV(flV-(bVW)))W'l = / drfW. (2.2) 
Jo. Jo. 

Integrating by parts, we obtain 

-^dT[(flV-(bVW))(flV-(bVW))] = j^dxfW, (2.3) 

where we have used the fact that W and its normal derivative vanish on 90. 
Importantly, we need to check that this problem is well posed. For this, we denote by 
L^(Q) (resp. L^(n)) the (class of) square integrable functions on O, with values in C 
(resp. C^). We also introduce the Hilbert space 

H^{n) = {v^L^{n.) : Vve\}{c^),^veL^{c^),v^ve\?{c^)} (2.4) 

If we further assume that W = and n - VW = on 30, the above space is denoted H3(0). 
The existence and uniqueness of the solution is then straightforwardly ensured by the 
Lax-Milgram theorem applied for / G L^(0), and {u,v) G Hq (O) x Hq (O), provided that 
there exist mi, m2, Mi, M2 > such that Mi >a>mi and M2 > K*hK > m2 for all k G C^. 



2.2 Discretisation of the variational equation 

To discretise this variational problem, the basic idea is to approximate W with a linear 
combination of simple functions (pj{x,y) called form functions and spanning a functional 
space of finite dimension (i.e. a Galerkin space with a finite number of degrees of free- 
dom). Then, this development is inserted into the weak form of the fourth-order partial 
differential equation thus generating a finite system of equations. 

It is common to consider a triangular mesh with a reference cell {x>0 ,y>0 : x+y < 
1} the discussion we assume a square mesh. We then consider e.g. a basis of second 
order polynomials {1, ,x,y , x^ ,y^ ,xy} which spans a space of dimension 6. High-order 
polynomial basis can also be used. The function W{x,y) can be written as 



N 

W = X:/3,</),(x,y) 



(2.5) 
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where /3, are values of u at the node points of the triangular mesh. The form functions 
satisfy the condition (pi(pj = Sij {S is the Kronecker symbol). 

Finally using the development of W (|2.5|), we obtain the system of n equations with n 
unknowns where we replaced W by the set of form functions (pi 

-L(^l^dT[{aV-{hV(Pj)){aV-{hVcPi))]^ = j^dTfcPi = G ,i = l,...,N , (2.6) 

and this can be written as a linear algebraic system K2A = Q which can be solved itera- 
tively. 



2.3 Perfecly Matched Layers for flexural waves 



R: 



R 



Incident flexural 
plane wave 




Arbitrary-shaped 
obstacle 



R: 



R: 



R. 



R: 



Figure 1: Schematic diagram showing the problem to be modelled: The plane elastic wave is incident from left 
to right in the presence of a diffracting obstacle of arbitrary shape. The different domains forming the PMLs 
regions are denoted respectively Ri, R2 and K3. 



Perfectly Matched Layers (PMLs) are equivalent to a geometrical transformation. In 
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Cartesian coordinates they can be expressed as 

f'X f'Xj f'Z 

Xs{x)= dx'sx{x') ys{y)= dy'sy{y') and Zs{z)= dz'sz{z'), (2.7) 
JO Jo Jo 

whose Jacobian is given by 

J„, = #^^ = Diag(iii) (2.8) 

where Diag represents a diagonal matrix. 

We consider now the transformation of the biharmonic equation under an arbitrary 
change of coordinates (say from {x,y,z) to (xs,t/s,Zs)). To do so, we have to consider the 
way the two equations in the following system transform: 

I V.(-r^Vy)+A-i^4^=0 ^^'^^ 

where ^ is an inhomogeneous anisotropic 2D tensor and A is an inhomogeneous coeffi- 
cient of the material (C^ = E^^''^ and \ = p^^^^). 

For example, the first equation of this system is similar to Helmholtz's equation except 
that here we have two functions V and W instead of a unique function. Though we can 
write it in an integral form in the domain of validity let's say O 

/ dTV.(-r^VW)+ / dT\-^V = Q (2.10) 
Jo. = Jci 

with dr the infinitesimal element of volume. 

By multiplying this equation by a test function \p and integrating by parts, we obtain 

-/ dTVtp-{-r^VW)+ f dTA-^tpV = (2.11) 
Jn = Jo. 

Now, by operating the previously described coordinate transform characterized by 
its jacobian ]x^x = Jxj, that we will will note simply / to lighten the notations, (|2.11l) trans- 
forms as 

where we have used the generic way the gradient and the infinitesimal volume trans- 
form: V = J* V and dr = dr' • Finally, using the useful expression of the scalar prod- 
uct A-B = A*B where A* is the transpose of A we can write this equation in the more 
adapted form 
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Figure 2: (a) Decay of the plane wave when it penetrates in the domain denoted by Ri or R2 in Figure [T] we 
represent here the amplitude of the plane wave in the direction i/ = or x = 0. (b) Propagation of the plane 
wave along the direction x = y when it lies in the central domain (left of the red vertical line) and when it lies 
in the PMLs domain R3 showing clearly that its amplitude vanishes rapidly in R3. 



Hence, we can see from this last expression, that transforming the biharmonic equa- 
tion from the initial system of coordinates to the new one is equivalent to replacing the 
initial parameters and A by the deduced parameters and A' 
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A'-i = A-Vdet(J,,,) 

We can obtain the same expressions by considering the second equation in (I2.9D . 

Now, consider that we have an isotropic and homogeneous media, the use of PMLs 
introduce the complex matrix Tp^^ which represents the inverse of the metric tensor. It 
can be expressed in terms of the parameters of the Jacobian of the PML transformation 

Sx, Sz and 

T,l,, = Diag(^,^,^) (2.15) 

^x °y ■^z 

We have 



, _1 , (^^_[ 1 In the domain of study 
s,-Sy-i,s,.^xj-| in the regional (See FigureUD 

For the regions R2, we have just to replace x by y and 1 by 2 in equation (|2.16|) , and 
for the regions J?3, we have 

s, = l, sy{y) = l-i^ andsx(x) = l-/-^ (2.17) 
^ eve eve 

Figure |2] clearly confirms the mechanism of PMLs based on the transformation of bi- 
harmonic equation from the initial system {x,y) to ixs,ys) representing absorbing layers 
in such a way incident waves are perfectly transmitted (without any reflexions) indepen- 
dently from their frequency or incidence angles. 

To our knowledge, this is the first time PMLs were applied to the fouth order biharmonic 
problem. 

For exapmle. Figure |2] (b) shows the decay of an incident plane wave in the PML zone 
(rectangle R3 of Figured]): the amplitude of the wave can be considered zero after only a 
short propagation distance in the absorber medium. 



3 Elliptic invisibility cloaks and the mirage effect for biharmonic 
problems 

3.1 Elliptic invisibility cloak for flexural waves 

This may be deduced from the circular case by scaling the Cartesian coordinates. The 
global transformation is a mapping of a holey elliptic domain (the inner and outer bound- 
aries are concentric ellipses with the same eccentricity) on a simply connected elliptic do- 
main bounded by the outer ellipse of the cloak. The detailed computation of the equiva- 
lent material properties is given by the following sequence of transformations. The start- 
ing point is the ellipse that bounds the exterior limit of our cloak with its principal axes 
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Figure 3: Right: Real part of the displacement field W distribution in the vicinity of the cloaked elliptic rigid 
clamped obstacle. The source is located at point (0.5,0.5) and its wavelength is Aq =27r//3Q = 0.28, which is 
of the same order as the inner and outer radii (major and minor) of the cloak: fli = 0.2, fl2 = 0.3 and fei = 0.4, 
^2 = 0.6. Left: Real part of the displacement field W scattered by a rigid clamped obstacle of minor radius 0.2 
m and major radius 0.3 m for an incoming cylindrical wave of the same frequency. 



chosen conveniently parallel to the coordinates axes. The first step is aimed at restoring 
the previous situation namely a circular cloak. For this purpose, the plane is scaled by a 
factor Sy along the (arbitrary chosen) -axis so that the initial ellipse becomes a circle (in 
the scaled coordinates) defined by the transformation y = Syi/s (and simply x = Xs) char- 
acterized by the Jacobian matrix ]xx, = Diag(l/Si/,1). The next three transformations are 
then the ones used to build the circular cloak: transformation to cylindrical coordinates, 
radial contraction (the active part), and transformation to rectangular coordinates. It is 
important to note that it is the scaled variable x/s that is involved in these various oper- 
ations. The last step to be performed is an inverse scaling along y-axis: y' = {1/Sy)y'g to 
recuperate the initial elliptical shape of the cloak and an identity for the transformation 
of the outside of the cloak. At the end, a cloak is obtained whose inner and outer bound- 
aries are ellipses with Ri = x-axis of the hole, Sy, Ri = y-axis of the hole, R2 = a:-axis of 
the external boundary, Sy, R2 = y-axis of the external boundary. The total Jacobian of this 
sequence of transformations is 

Jxx' — J xxg] Xsr] rr'] r' x'g] x'^x' (3.1) 

The inverse of the matrix is given explicitly by ||5l 

T-i = Diag(l,Sy,l)R(0ODiag(^^,^,l)R(-0,)Diag(l,^,l)R(0O 
Diag(a,4,l)R(-0ODiag(l,Sy,l)^, 
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with r', = ^Jx'^ + {y'/Syf, Os = Q's = 2arctan( (y'/sy ) (x' +4) ) and r = (r^ - Ri ) / a. 

Note the angles and distances are computed in the scaled coordinate systems {xs,ys) 
where the ellipses are mapped to circles. 

Figure|3](b) shows the calculated displacement field distribution when an elastic point 
source vibrating in the out-of-plane direction is placed near the optimized cloak. It is 
shown that in the near region of the cloak, the displacement field stays almost unper- 
turbed, thus a so called perfect elliptic cloak is obtained. In order to see how much 
improvement has been achieved through the presence of the elliptic coat, snapshots of 
the displacement W in the presence of a clamped obstacle ( W = and 3 W/3r = where r 
describes the boundary of the ellipse) are shown in (a). 

3.2 Transformation elastodynamics: Elliptic rotator and mirage effect 




Figure 4: A plane bending wave is incident from left to right onto a coating consisting of two concentric ellipses 
of the same dimensions as in Figure [3] The Young modulus and density of the shell are deduced from the 
transformation given in Equation l|3.3p . We can see that the wave fronts are rotated with an angle f?o = 7r/4 
radians. 

Based on the coordinate transformation methods [18J, we propose here a two-dimensional 
transformation-media coat that can rotate fields but remains itself invisible. The Young 
modulus tensor and density of such a "field rotator" for flexural waves can be deduced 
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from appropriate coordinates transforms. 

The difference with invisibility (considered in the previous section and where we map 
a point on a circle of radius Ri ) is that here, we map a circle (of radius Ki or R2) onto 
itself. So let's consider the following transformation: 



with oci = 6q/ (Ri — R2) and ^i = 9q/ (R2 — Ri). Oq is the angle of rotation of the wave fronts. 
The Jacobian of this transformation is given by 



This expression may be injected in Equation (13. II) instead of the invisibility Jacobian Jyyi to 
obtain the matrix T which gives us the physical parameters of the coat (Young's modulus 
and density). In addition, inside the inner circle of radius Ri, we have the additional 
transform (r' = r, 9' = 6+0q and z' = z) whose Jacobian is the 3x3 identity matrix and 
doesn't contribute thus to the materials specifications inside the cloak. Finally, the FEM 
computations using these expression show in Figure|4]the displacement filed in the vicin- 
ity of an ideal field rotator designed by transformation-media theory. 

We consider also another application of coordinates transformation and which has 
been verified for others types of waves (acoustic and elastodynamics). 

In Figure |5] the vertical displacement W outside the cloak seems to originate from a 
location r^, which is slightly shifted with respect to the real position of the source (lo- 
cated inside the cloak itself). This effect is similar to the mirage effect already observed in 
electromagnetic cloaks [3J. The vertical displacement on the line x = y is shown for both 
a source located inside the coating and a source shifted laterally at r, in a homogeneous 
plate (i.e., without any cloak and inclusion). The respective positions are given in the 
figure caption. 

In conclusion, we have studied in this paper some properties of the biharmonic equa- 
tion governing the propagation of flexural waves in thin elastic plates, we have shown 
the way to generalize some concepts like cloaking or the rotator and the mirage effect 
to the domain of elasticity in the special case of elliptic geometries. Despite the fact that 
general elastic equations are not invariant under geometrical changes, we have demon- 
strated that cloaking through coordinates changes is applicable to some special configu- 
rations (thin plates). Moreover Perfectly Matched Layers have been performed using the 
finite elements technique for the biharmonic problems. 




(3.3) 




(3.4) 
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Figure 5: (a) Real part of the displacement field W when the source is located inside the coating at the point 
(0.25,0.25). The field seems to be emitted by a shifted source located at point (0.425,0.325) (b); (c) Real 
part of the displacement W along the line x = y for a source located inside the coating at point (0.25,0.25) as 
shown in Figure [5] (a) (solid line) and for a source located in the homogeneous plate at point (0.425,0.325) 
(dashed-starred line) deduced from the geometrical transform given in l|3.2p . 



Acknowledgments 



The authors would like to thank insightful discussions with Prof. A.B. Movchan and 
N.V. Movchan at Liverpool University, Dr. C.G. Poulton at the University of Technology 



14 



of Sydney (UTS) and Prof. R.C. McPhedran at the University of Sydney. 
References 

[1] J.B. Pendry, D. Schurig, and D.R. Smith, Science 312, 1780 (2006). 

[2] D. Schurig, J.J. Mock, B.J. Justice, S.A. Cummer, J.B. Pendry, A.E Starr, D.R. Smith, Science 
314, 977 (2006). 

[3] F. Zolla, S. Guenneau, A. Nicolet, and J.B. Pendry, Opt. Lett. 32, 1069-1071 (2007). 

[4] W. Cai, U.K. Chettiar, A.V. Kildiev and V.M. Shalaev, Nature 1, 224-227 (2007). 

[5] A. Nicolet, F. Zolla and S. Guenneau, IEEE Trans. MAg. 44, 1150-1153 (2008) 

[6] S.A. Cummer and D. Schurig, New J. Phys. 9, 45 (2007). 

[7] D. Torrent and J. Sanchez-Dehesa, New J. Phys. 10, 063015 (2008). 

[8] M. Farhat, S. Enoch, S. Guenneau and A.B. Movchan, Phys. Rev. Lett. 101, 134501 (2008). 

[9] H. Chen and C. T. Chan, Appl. Phys. Lett. 91, 183518 (2007). 
[10] G.W. Milton, M. Briane, and J.R. Willis, New J. Phys. 8, 248 (2006). 
[11] M. Brun, S. Guenneau and A.B. Movchan, Appl. Phys. Lett. 94 061903 (2009). 
[12] H. Chen, B.I. Wu, B. Zhang, and J.A. Kong, Phys. Rev Lett. 99, 063903 (2007). 
[13] U. Leonhardt and T G. Philbin, New J. Phys. 8, 247 (2006). 

[14] A. Nicolet, J.E Remade, B. Meys, A. Genon and W. Legros, Appl. Phys. 75, 6036-6038 (1994). 
[15] F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, S. Guenneau and D. Felbacq, Foundations of 

photonic crystal fibres (Imperial College Press, London, 2005). 
[16] E.G. Post, Formal Structure of Electromagnetics; General Covariance and Electromagnetics 

(Interscience, 1962). 
[17] A.J. Ward and J.B. Pendry J. Mod. Opt. 43, 773-793 (1996). 
[18] U. Leonhardt, Science 312 1777-1780 (2006). 

[19] N.A. Nicorovici, R.C. McPhedran and G.W. Milton, Phys. Rev B 49, 8479-8482 (1994). 
[20] J.B. Pendry Phys. Rev Lett. 86, 3966-3969 (2000). 
[21] A. Alu and N. Engheta, Phys. Rev E 95 016623 (2005). 

[22] B. Edwards, A. Alu, M. G. Silveirinha, and N. Engheta, Phys. Rev Lett. 103 153901 (2009). 
[23] M. Farhat, S. Guenneau, A.B. Movchan and S. Enoch, Opt. Express 16, 5656-5661 (2008) 
[24] S.A. Cummer, B.I. Popa, D. Schurig, D.R. Smith, J. Pendry, M. Rahm, and A. Starr, Phys. Rev 

Lett. 100, 024301 (2008). 
[25] M. Farhat, S. Enoch, S. Guenneau and A.B. Movchan, Phys. Rev B 79 033102 (2009). 
[26] S. Timoshenko, Theory of plates and shells (McGraw-Hill, New York, 1940). 
[27] K.F. Graff, Wave motion in elastic solids (Dover, New York, 1975). 
[28] L.D. Landau, and E.M. Lifschitz, Elasticity theory (Pergamon Press, 1954). 



